#ifndef CLASS_GF1IO
#define CLASS_GF1IO

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <sys/stat.h>  // for file sizeof

#include "Timer.h"


class GF1IO {
	public:
	
	//P :: VS Gadget/Arepo Particle Data Container
	class P { 
	public:
		float  Pos[3];
		float  Vel[3];
		float  Mass;
		uint 	 Type;
		float  pot;
		
		float  Rho, U, Temp, Ne;
	};
		
	class Header {
		public:
		Header():
		num_files(1)
		{}
		
		int      npart[6];
		double   mass[6];
		double   time;
		double   redshift;
		int      flag_sfr;
		int      flag_feedback;
		unsigned int      npartTotal[6];
		int      flag_cooling;
		int      num_files;
		double   BoxSize;
		double   Omega0;
		double   OmegaLambda;
		double   HubbleParam; 
		
		int flag_stellarage;		/*!< flags whether the file contains formation times of star particles */
		int flag_metals;		/*!< flags whether the file contains metallicity values for gas and star
				   particles */
		
		unsigned int npartTotalHighWord[6];	/*!< High word of the total number of particles of each type */
		int flag_entropy_instead_u;	/*!< flags that IC-file contains entropy instead of u */
		int flag_doubleprecision;	/*!< flags that snapshot contains double-precision instead of single precision */

		int flag_ic_info;             /*!< flag to inform whether IC files are generated with ordinary Zeldovich approximation,
                                     or whether they ocontains 2nd order lagrangian perturbation theory initial conditions.
                                     For snapshots files, the value informs whether the simulation was evolved from
                                     Zeldoch or 2lpt ICs. Encoding is as follows:
                                        FLAG_ZELDOVICH_ICS     (1)   - IC file based on Zeldovich
                                        FLAG_SECOND_ORDER_ICS  (2)   - Special IC-file containing 2lpt masses
                                        FLAG_EVOLVED_ZELDOVICH (3)   - snapshot evolved from Zeldovich ICs
                                        FLAG_EVOLVED_2LPT      (4)   - snapshot evolved from 2lpt ICs
                                        FLAG_NORMALICS_2LPT    (5)   - standard gadget file format with 2lpt ICs
                                     All other values, including 0 are interpreted as "don't know" for backwards compatability.
                                 */
		float lpt_scalingfactor;      /*!< scaling factor for 2lpt initial conditions */
		
		char fill[256- 6*4- 6*8- 2*8- 2*4- 6*4- 2*4 - 4*8 - 4*2 - 6*4 - 3*4 - 1*4];  /* fills to 256 Bytes */
	
	} header;
	
	int NumPart, Ngas;
	
	P *p[6];
	
	uint *Id;
	uint blocksize;
	uint pTypeNum;
	char *mem;
	char *filename, *filename_base;
	uint *ofs;
	
	class TN {
		public:
			int N[7];
			int &operator =(int i)   {  N[0]=i; return i;  }
			int &operator [](int i)  {  return N[i+1];     }
			operator int(){	return N[0]; }
	} N;
	
	GF1IO():
	Id(0),
	blocksize(4),
	pTypeNum(6),
	mem(0),
	filename(0),
	filename_base(0),
	ofs(0)
	{ *p=0; }
	
	~GF1IO(){
		if (*p) delete [] *p;
		if (mem) delete [] mem;
		if (Id) delete [] Id;
	}
	
	
	void free(){
		
		if (*p) delete [] *p;
		if (mem) delete [] mem;
		if (Id) delete [] Id;
		
		*p=0; mem=0; Id=0;
	}

	
	
	double Time, Redshift;
	double readtime;
	double readspeed;
	uint filesize;
	int N_withVarMassess;
	
	int getFileSize(const char* fullfilename){
		struct stat results;
		stat(fullfilename, &results);
		return results.st_size;
	}

	
	char* readraw(const char* filename) {
		
		filesize = getFileSize(filename);
		char *mem = new char[filesize];
		
		Timer tr; 
		FILE *fs = fopen(filename, "r");
		if (fs==NULL) {printf("\033[22;31merror: can't read the file\033[22;35m '%s'\033[22;30m\n", filename); return NULL;}
		
		tr.set();
		fread(mem, 1, filesize, fs); 
		tr.get();
		
		readspeed = filesize/tr.sin;
		
		fclose(fs);
		
		return mem;
		
	}

	char* getSplitFileName(const char *basefilename, uint ifile){
		char *splitfilename = new char[256];
		sprintf(splitfilename, "%s.%u", basefilename, ifile);
		return splitfilename;
	}
	
	void read(const char* filename){
		
		
		
		Timer tr; 
		
		//fnd.decompose(filename);
		
		if (this->filename) delete [] this->filename;
		this->filename = new char[256];
		sprintf(this->filename, "%s", filename);
		
		uint ifile=0;
		
		char *mem0;
		
		uint *ofs = new uint[pTypeNum];
		memset(ofs,0,pTypeNum*sizeof(uint));
		
		uint filesize(0);
		
		do {
			
		
		if (mem)	delete [] mem;
		tr.set();
		mem0 = mem = readraw(this->filename);
		
		tr.get();
		filesize+=this->filesize;
		
		//header
		mem+=blocksize;
		memcpy(&header, mem, sizeof(Header)); 
		mem+=sizeof(Header);
		mem+=blocksize;
		
		
		// local
		uint NumPart=0;
		uint Ngas= header.npart[0];
		for (uint i=1; i<pTypeNum; i++)
		NumPart+= header.npart[i];
		uint N=NumPart+Ngas;
		
		// global
		this->NumPart=0;
		this->Ngas= header.npartTotal[0];
		for (uint i=1; i<pTypeNum; i++)
		this->NumPart+= header.npartTotal[i];
		
		
		this->N=this->NumPart+this->Ngas;
		for (uint i=0; i<pTypeNum; i++)
		this->N[i] = header.npartTotal[i];
		
		
		
		// offsets
		if (ifile==0) {
			if (this->ofs) delete [] this->ofs;
			this->ofs = new uint[pTypeNum];
			memset(this->ofs,0,pTypeNum*sizeof(uint));
			for (uint k=0; k<pTypeNum-1; k++)
			this->ofs[k+1] = this->ofs[k] + header.npartTotal[k];
		}

		
		// p
		if (!ifile) {
			if (*p) delete [] *p;
			*p = new P[this->N];
		}
		
		
		
		for (uint i=1; i<pTypeNum; i++)
		p[i] = p[i-1] + this->N[i-1];
		



		// Pos
		mem+=blocksize;
		for (uint k=0; k<pTypeNum; k++)
		for (uint j=0; j<header.npart[k];j++) {
		memcpy((*p)[ this->ofs[k]+ofs[k]+j ].Pos, mem, 3*sizeof(float));  
		mem+=3*sizeof(float); 
		}
		mem+=blocksize;
		
		
		
		// Vel
		mem+=blocksize;
		for (uint k=0; k<pTypeNum; k++)
		for (uint j=0; j<header.npart[k];j++) {
		memcpy((*p)[ this->ofs[k]+ofs[k]+j ].Vel, mem, 3*sizeof(float));  
		mem+=3*sizeof(float); }
		mem+=blocksize;
		
		// Id
		if (!ifile) {
			if (Id) delete [] Id;
			Id = new uint[this->N];
		}
		
		// Ids
		mem+=blocksize;
		for (uint k=0; k<pTypeNum; k++) {
		memcpy(&Id[this->ofs[k]+ofs[k]], mem, header.npart[k]*sizeof(uint));  
		mem+=header.npart[k]*sizeof(uint);
		}
		mem+=blocksize;

		// Mass
		N_withVarMassess=0; // lets check whether the var masses block is present
		for (uint k=0; k<pTypeNum; k++)
		if (header.mass[k]==0)
		N_withVarMassess+=header.npart[k];
		
		
		
		if (N_withVarMassess>0) mem+=blocksize; // do it if the var masses block is present
		for (uint k=0; k<pTypeNum; k++) 
		for (uint i=0; i<header.npart[k]; i++) {
			if (header.mass[k]==0) { // we read the block
				memcpy(&((*p)[this->ofs[k]+ofs[k]+i].Mass), mem, 1*sizeof(float));  
				mem+=1*sizeof(float); 
			}
			else // in this case we use common masses
				(*p)[this->ofs[k]+ofs[k]+i].Mass = header.mass[k];
			
		}
		if (N_withVarMassess>0) mem+=blocksize;

		
		// Gas
		if(header.npart[0]>0) {
			
			// Internal Energy
			mem+=blocksize;
			for (uint j=0; j<header.npart[0]; j++) {
			memcpy(&(*p)[*ofs+j].U, mem, 1*sizeof(float));  
			mem+=1*sizeof(float);	}
			mem+=blocksize;

			// Density
			mem+=blocksize;
			for (uint j=0; j<header.npart[0]; j++) {
			memcpy(&(*p)[*ofs+j].Rho, mem, 1*sizeof(float));  
			mem+=1*sizeof(float);  }
			mem+=blocksize;
		
			
			if(header.flag_cooling) { // this is for the cooling in the user guid it says it is currently unused
				mem+=blocksize;
				for (uint j=0; j<header.npart[0]; j++) {
					memcpy(&(*p)[*ofs+j].Ne, mem, sizeof(float));
					mem+=sizeof(float);
				}
				mem+=blocksize;
			} else 
				for (uint j=0; j<header.npart[0]; j++)
				(*p)[*ofs+j].Ne= 1.0;
		
		}
		
		
		for (uint k=0; k<pTypeNum; k++) 
		ofs[k]+=header.npart[k];
		
		mem = mem0;	
		
		if (header.num_files>1) { 
			
			ifile++; 
			
			if (this->filename) delete [] this->filename;
			this->filename = getSplitFileName(filename, ifile);
			
		} else break;
		
		} while (ifile<header.num_files);
 
		Time= header.time; // read snapshot time (or exspansion factor) 
		Redshift= header.time; // read the snapshot redshift 
		
		readspeed=filesize/tr.sin_tot;
		readtime=tr.sin_tot;
		this->filesize=filesize;
		delete [] ofs;


	}
	
	void write(const char *filename){
		
		printf("write::#00\n");
		
		char *mem0(0);
		
		uint *ofs = new uint[pTypeNum];
		memset(ofs,0,pTypeNum*sizeof(uint));
		
		// I need offset
		if (this->ofs) delete [] this->ofs;
		this->ofs = new uint[pTypeNum];
		*(this->ofs)=0;
		for (uint k=0; k<pTypeNum-1; k++)
		this->ofs[k+1] = this->ofs[k] + header.npartTotal[k];
		
		printf("write::#01\n");
		
		printf("header.num_files=%i\n", header.num_files);
		
		if (header.num_files==0) header.num_files = 1;
			
 		for (uint ifile=0; ifile<header.num_files; ifile++){
		
		
			for (uint k=0; k<pTypeNum; k++) 
			if (header.npart[k]) {
				header.npart[k]=header.npartTotal[k]/header.num_files;
				if (ifile==header.num_files-1) header.npart[k]+=header.npartTotal[k]%header.num_files;
			} else header.npart[k]=0;
		
			
			uint n_tot(0);
			for (uint k=0; k<pTypeNum; k++) 
			if (header.npart[k]) n_tot+=header.npart[k];
		
			
			
			uint n_totvm(0);
			for (uint k=0; k<pTypeNum; k++) 
			if (header.npart[k] && !header.mass[k]) n_totvm+=header.npart[k];
			
			
			// memory size defenition
			filesize = sizeof(Header) 
			+  n_tot* sizeof( *Id ) 
			+	n_tot*3*( sizeof( *(**p).Pos ) +  sizeof( *(**p).Vel ) ) 
			+	(1+1 + 1+1 + 1+1 + 1+1)*blocksize;
			
			for (uint k=0; k<pTypeNum; k++)
			if (header.mass[k]==0) filesize+= header.npart[k]*sizeof((**p).Mass);
			
			if (n_totvm) filesize+=2*blocksize;
			if (header.npart[0]) filesize+=4*blocksize + n_tot*(sizeof((**p).U) + sizeof((**p).Rho));
			if (header.flag_cooling) filesize+=2*blocksize + n_tot*sizeof((**p).Ne);
			
			if (mem0) delete [] mem0;
			mem0 = mem = new char[filesize];

			
			
			*(uint*)mem = sizeof(Header);
			mem+=blocksize;
			memcpy(mem, &header, sizeof(Header));
			mem+=sizeof(Header);
			*(uint*)mem = sizeof(Header);
			mem+=blocksize;
			
			
			// Pos
			*(uint*)mem = 3*n_tot*sizeof(float);
			mem+=blocksize;
			for (uint k=0; k<pTypeNum; k++)
			for (uint j=0; j<header.npart[k];j++) {
				memcpy(mem, (*p)[ this->ofs[k] + ofs[k] + j ].Pos, 3*sizeof(float));  
				mem+=3*sizeof(float); 
			}
			*(uint*)mem = 3*n_tot*sizeof(float);
			mem+=blocksize;
		
			
			
			// Vel
			*(uint*)mem = 3*n_tot*sizeof(float);
			mem+=blocksize;
			for (uint k=0; k<pTypeNum; k++)
			for (uint j=0; j<header.npart[k];j++) {
				memcpy(mem, (*p)[ this->ofs[k] + ofs[k] + j ].Vel, 3*sizeof(float));  
				mem+=3*sizeof(float); 
			}
			*(uint*)mem = 3*n_tot*sizeof(float);
			mem+=blocksize;
			
			// Ids
			*(uint*)mem = n_tot*sizeof(float);
			mem+=blocksize;
			for (uint k=0; k<pTypeNum; k++) {
				memcpy(mem, &Id[this->ofs[k]+ofs[k]], header.npart[k]*sizeof(uint));  
				mem+=header.npart[k]*sizeof(uint);
			}
			*(uint*)mem = n_tot*sizeof(float);
			mem+=blocksize;
			
			// Var Masses
			if (n_totvm) {
				*(uint*)mem = n_totvm*sizeof(float);
				mem+=blocksize;
				for (uint k=0; k<pTypeNum; k++) 
				for (uint i=0; i<header.npart[k]; i++)
				if (header.mass[k]==0) { 
					memcpy(mem, &(*p)[this->ofs[k]+ofs[k]+i].Mass, sizeof(float));  
					mem+=sizeof(float); 
				}
				*(uint*)mem = n_totvm*sizeof(float);
				mem+=blocksize;
			}
			
			// Gas
			if (header.npart[0]) {
				
				
				// Internal Energy
				*(uint*)mem = header.npart[0]*sizeof(float);
				mem+=blocksize;
				for (uint j=0; j<header.npart[0]; j++) {
				memcpy(mem, &(*p)[*ofs+j].U, sizeof(float));  
				mem+=sizeof(float);	}
				*(uint*)mem = header.npart[0]*sizeof(float);
				mem+=blocksize;
				
				
				// Density
				*(uint*)mem = header.npart[0]*sizeof(float);
				mem+=blocksize;
				for (uint j=0; j<header.npart[0]; j++) {
				memcpy(mem, &(*p)[*ofs+j].Rho, sizeof(float));  
				mem+=sizeof(float);	}
				*(uint*)mem = header.npart[0]*sizeof(float);
				mem+=blocksize;
				
				// Cooling
				if (header.flag_cooling) { // this is for the cooling in the user guid it says it is currently unused
					*(uint*)mem = header.npart[0]*sizeof(float);
					mem+=blocksize;
					for (uint j=0; j<header.npart[0]; j++) {
					memcpy(mem, &(*p)[*ofs+j].Ne, sizeof(float));
					mem+=sizeof(float);	}
					*(uint*)mem = header.npart[0]*sizeof(float);
					mem+=blocksize;
				} 
				
			}
			
			printf("file[%u]:\n", ifile);
			char ifilename[256];
			header.num_files>1 ?
			sprintf(ifilename, "%s.%u", filename, ifile) :
			sprintf(ifilename, "%s", filename);
			
			
			FILE *pFile = fopen (ifilename , "wb" );
			fwrite(mem0, sizeof(char), filesize,  pFile);
			fclose(pFile);
			
			
			for (uint k=0; k<pTypeNum; k++) 
			ofs[k]+=header.npart[k];
			
			printf("file[%u] is passed.\n\n", ifile);
		}
		
		/*for (uint k=0; k<pTypeNum; k++)
		printf("ofs[%u]: %u\n", k, ofs[k]);*/
		mem = mem0;
		delete [] ofs;
	}
	
	
	
	
	void print_P(uint nprint){
		
		for (uint k=0, ofs=0; k<pTypeNum; ofs+=header.npartTotal[k], k++) {
			printf("%u: \033[22;37m%u, %u\033[22;30m\n", k, ofs, header.npartTotal[k]);
			
			if (header.npartTotal[k]) 
				for (uint j=0; j<=nprint; j++){
					uint ik = (header.npartTotal[k]-1)*double(j)/nprint; // becouse double precision problem
					//uint ik = j;
					uint i = ofs+ik;
					
					printf("\033[22;34mP\033[22;30m[%7u]:  \033[22;35m%10u\033[22;30m  % .7e", i, Id[i], (*p)[i].Mass);
					printf("\t% .7e % .7e % .7e", 	(*p)[i].Pos[0], (*p)[i].Pos[1], (*p)[i].Pos[2]);
					printf("\t% .7e % .7e % .7e", 	(*p)[i].Vel[0], (*p)[i].Vel[1], (*p)[i].Vel[2]);
					if (k==0 && header.npartTotal[k]>0) 
					printf("\t% .7e % .7e % .7e", 	(*p)[i].U, (*p)[i].Rho, (*p)[i].Ne);
					printf("\t\033[22;36m %u \033[22;30m", (*p)[i].Type);
					printf("\n");
					
				}
		}
		
	}
	
	void print(){
		
		const double MB = 1<<20;
		
		printf("\n");
		printf("\033[22;34mFilename\033[22;30m: \033[22;32m'%s'\033[22;30m\n", filename);
		
		printf("\033[22;34mFilesize\033[22;30m: %.3f MB \033[22;37m( %u )\033[22;30m\n", filesize/MB, filesize);
		printf("\033[22;34mReed Speed\033[22;30m: %.2f MB/sec\n", readspeed/MB);
		printf("\033[22;34mReed Time\033[22;30m: %.2f sec\n", readtime);
		printf("\033[22;34mNumFilesPerSS\033[22;30m: %u\n\n", header.num_files);
		
		printf("\033[22;34mBoxSize\033[22;30m: %.15g\n", header.BoxSize);
		printf("\033[22;34mOmega0\033[22;30m: %.15g\n", header.Omega0);
		printf("\033[22;34mOmegaLambda\033[22;30m: %.15g\n", header.OmegaLambda);
		printf("\033[22;34mHubbleParam\033[22;30m: %.15g\n", header.HubbleParam);
		printf("\033[22;34mRedshift\033[22;30m: %.15g\n", header.redshift);
		printf("\033[22;34mTime\033[22;30m: %.15g\n", header.time);

		int flag_stellarage;		/*!< flags whether the file contains formation times of star particles */
		int flag_metals;		/*!< flags whether the file contains metallicity values for gas and star
				   particles */
		
		unsigned int npartTotalHighWord[6];	/*!< High word of the total number of particles of each type */
		int flag_entropy_instead_u;	/*!< flags that IC-file contains entropy instead of u */
		int flag_doubleprecision;	/*!< flags that snapshot contains double-precision instead of single precision */

		int flag_ic_info;             /*!< flag to inform whether IC files are generated with ordinary Zeldovich approximation,
                                     or whether they ocontains 2nd order lagrangian perturbation theory initial conditions.
                                     For snapshots files, the value informs whether the simulation was evolved from
                                     Zeldoch or 2lpt ICs. Encoding is as follows:
                                        FLAG_ZELDOVICH_ICS     (1)   - IC file based on Zeldovich
                                        FLAG_SECOND_ORDER_ICS  (2)   - Special IC-file containing 2lpt masses
                                        FLAG_EVOLVED_ZELDOVICH (3)   - snapshot evolved from Zeldovich ICs
                                        FLAG_EVOLVED_2LPT      (4)   - snapshot evolved from 2lpt ICs
                                        FLAG_NORMALICS_2LPT    (5)   - standard gadget file format with 2lpt ICs
                                     All other values, including 0 are interpreted as "don't know" for backwards compatability.
                                 */
		float lpt_scalingfactor;      /*!< scaling factor for 2lpt initial conditions */
		
		
		printf("\033[22;34mflag_stellarage\033[22;30m: %i\n", header.flag_stellarage);
		printf("\033[22;34mflag_metals\033[22;30m: %i\n", header.flag_metals);
		printf("\033[22;34mflag_entropy_instead_u\033[22;30m: %i\n", header.flag_entropy_instead_u);
		printf("\033[22;34mflag_doubleprecision\033[22;30m: %i\n", header.flag_doubleprecision);
		printf("\033[22;34mflag_ic_info\033[22;30m: %i\n", header.flag_ic_info);
		printf("\033[22;34mlpt_scalingfactor\033[22;30m: %g\n", header.lpt_scalingfactor);
		
		
		
		printf("\033[22;34mNumPart\033[22;30m: %u\n", NumPart);
		printf("\033[22;34mNgas\033[22;30m: %u\n", Ngas);
		printf("\033[22;34mN\033[22;30m: \033[22;35m%u\033[22;30m\n\n", Ngas+NumPart);
		
		for (uint i=0; i<pTypeNum; i++)
			header.npart[i] ?
			printf("\033[22;34mnpart\033[22;30m[%u] = \033[22;33m%u\033[22;30m\n", i, header.npart[i]) :
			printf("\033[22;37mnpart[%u] = %u\033[22;30m\n", i, header.npart[i]);
			printf("\n");

		for (uint i=0; i<pTypeNum; i++)
			header.npartTotalHighWord[i] ?
			printf("\033[22;34mnpartTotalHighWord\033[22;30m[%u] = \033[22;33m%u\033[22;30m\n", i, header.npartTotalHighWord[i]) :
			printf("\033[22;37mnpartTotalHighWord[%u] = %u\033[22;30m\n", i, header.npartTotalHighWord[i]);
			printf("\n");


		
		for (uint i=0; i<pTypeNum; i++)
			header.mass[i] ?
			printf("\033[22;34mmass\033[22;30m[%u] = %.15g\n", i, header.mass[i]) :
			printf("\033[22;37mmass[%u] = %g\033[22;30m\n", i, header.mass[i]);
			printf("\n");
		
		for (uint i=0; i<pTypeNum; i++)
			header.npartTotal[i] ?
 			printf("\033[22;34mnpartTotal\033[22;30m[%u] = \033[22;33m%u\033[22;30m\n", i, header.npartTotal[i]) :
			printf("\033[22;37mnpartTotal[%u] = %u\033[22;30m\n", i, header.npartTotal[i]);
			printf("\n");
		
	}
	
	
	
	
};

#endif


// MODIFICATIONS

/* 1. In version .02 was memory consumption problem after deconstruction of GSIO object.
 * The reason was after deletion on pointers to arrays erased and that is way you cant 
 * to free memeory from prevous using like this: if (p) delete [] p
 * 
 * 2. there was also some problem with built-in baseName Finder, so i replased it with FND,
*/


  











